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Abstract 

The solar convection zone exhibits a strong level of differential rotation, whereby the 
rotation period of the polar regions is about 25-30% longer than the equatorial regions. The 
Coriolis force associated with these zonal flows perpetually "pumps" the convection zone 
fluid, and maintains a quasi-steady circulation, poleward near the surface. What is the 
influence of this meridional circulation on the underlying radiative zone, and in particular, 
does it provide a significant source of mixing between the two regions? In Paper I, we began 
to study this question by assuming a fixed meridional flow pattern in the convection zone and 
calculating its penetration depth into the radiative zone. We found that the amount of mixing 
caused depends very sensitively on the assumed flow structure near the radiative-convective 
interface. We continue this study here by including a simple model for the convection zone 
"pump", and calculating in a self-consistent manner the meridional flows generated in the 
whole Sun. We find that the global circulation timescale depends in a crucial way on two 
factors: the overall stratification of the radiative zone as measured by the Rossby number 
times the square root of the Prandtl number, and, for weakly stratified systems, the presence 
or absence of stresses within the radiative zone capable of breaking the Taylor-Proudman 
constraint. We conclude by discussing the consequences of our findings for the solar interior 
and argue that a potentially important mechanism for mixing in Main Sequence stars has 
so far been neglected. 



1. Introduction 

Various related mechanisms are thought to contribute to the generation and maintenance 
of large-scale meridional flows in the solar convection zone. The effect of rotation on turbulent 
convection induces a relatively strong anisotropy in the Reynolds stresses (Kippenhahn 1963), 
in particular near the base of the convection zone where the convective turnover time is of the 
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order of the solar rotation period. The divergence of these anisotropic stresses can directly 
drive large-scale meridional flows (see Rudiger, 1989, for a discussion of this effect). It also 
drives large-scale zonal flows (more commonly referred to as differential rotation) which then 
induce meridional forcing through the biais of the Coriolis force, a mechanism referred to as 
"gyroscopic pumping" (Mclntyre 2007). Indeed, the polar regions of the convection zone are 
observed to be more slowly rotating than the bulk of the Sun (Schou et al. 1998), so that 
the associated Coriolis force in these regions drives fluid towards the polar axis. Meanwhile, 
equatorial regions are rotating more rapidly than the average, and are therefore subject 
to a Coriolis force pushing fluid away from the polar axis. The most likely flow pattern 
resulting from the combination of these forces is one with an equatorial upwelling, a surface 
poleward flow and a deep return flow. This pattern is indeed observed near the solar surface: 
poleward surface and sub-surface flows with velocities up to a few tens of meters per second 
have been observed by measurements of photospheric line-shifts (Labonte & Howard 1982) 
and by time-distance helioseismology (Giles et al. 1997) respectively. 

The amplitude and spatial distribution of these meridional flows deeper in the convection 
zone remains essentially unknown, as the sensitivity of helioseismic methods rapidly drops 
below the surface. As a result, the question of whether some of the pumped mass flux actually 
penetrates into the underlying radiative zone is still open, despite its obvious importance for 
mixing of chemical species (Pinsonneault, 1997; Elliott & Gough, 1999), and its presumed 
role in the dynamical balance of the solar interior (Gough & Mclntyre 1998, Mclntyre, 
2007, Garaud, 2007, Garaud & Garaud, 2008) and in some models of the solar dynamo (see 
Charbonneau, 2005 for a review). 

In Paper I (Garaud & Brummell, 2008), we began a systematic study of the penetration 
of meridional flows from the convection zone into the radiative zone by considering a related 
but easier question: assuming that the amplitude and geometry of meridional flows in the 
convection zone are both known, what is their influence on the underlying radiative zone? 
This simpler question enabled us to study the dynamics of the radiative zone only by assum- 
ing a flow profile at the radiative-convective interface (instead of having to include the more 
complex convection zone in the calculation). The overwhelming conclusion of that first study 
was that the degree to which flows penetrate into the stratified interior (in the model) is very 
sensitive to the interfacial conditions selected. Hence, great care must be taken when using 
a "radiative-zone-only" model to make definite predictions about interior flow amplitudes. 
In addition, that approach makes the implicit assumption that the dynamics of the radiative 
zone do not in return influence those of the convection zone, but the only way to verify this 
is to construct a model which includes both regions. This was the original purpose of the 
present study; as we shall see, the combined radiative-convective model we construct here 
provides insight into a much broader class of problems. 
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We therefore propose a simplified model of the Sun which includes both a "connective" 
region and a "radiative" region, where the convective region is forced in such a way as to 
promote gyroscopic pumping of meridional flows. We calculate the flow solution everywhere 
and characterize how it scales in terms of governing parameters (e.g. stellar rotation rate, 
stratification, diffusivities, etc.), focusing in particular on the flows which are entering the 
radiative zone. We begin with a simple Boussinesq Cartesian model (Section [2]), first in the 
unstratified limit (Section 12. 3j) and then in the more realistic case of a radiative-convective 
stratification (Section 12. 4p . Although the Cartesian results essentially illustrate most of the 
relevant physical phenomena, we confirm our analysis with numerical solutions of the full set 
of equations in a spherical geometry in Section [3j We then use this information in Section 
H] to discuss the effects of mixing by meridional flows both in the Sun and in other Main 
Sequence stars. 

2. A Cartesian model 
2.1. Model setup 

As in Paper I, we first study the problem in a Cartesian geometry. Since our primary aim 
is to understand the behavior of the meridional flows generated (e.g. scaling of the solutions) 
in terms of the governing parameters, this approach is sufficient and vastly simplifies the 
required algebra. In Section [3], we turn to numerical simulations to study the problem in a 
spherical geometry. 

In this Cartesian model section, distances are normalized to the solar radius R Q , and 
velocities to -R©f2© where Q Q is the mean solar angular velocity (the exact value is not 
particularly relevant here). The coordinate system is (x,y,z), where x should be thought of 
as the azimuthal coordinate 0, with x G [0, 2ir]\ y represents minus the co-latitude and spans 
the interval y G [0, tc] (the poles are at y — and y = ir while the equator is at y = vr/2). 
Finally the z-direction is the radial direction with z G [0, 1], and represents the direction of 
(minus) gravity so that z = is the interior, and z — 1 is the surface. 

In this framework, the system rotates with mean angular velocity fl = (0, 0, 1), thereby 
implicitly assuming that the rotation axis is everywhere aligned with gravity. This assump- 
tion induces another "geometric" error in the velocity estimates for the meridional flows, 
comparable with the error made in reducing the problem to a Cartesian analysis; it does not 
influence the predicted scalings (except in small equatorial regions which we ignore here). 

We divide the domain in two regions, by introducing the dimensionless constant h to 
represent the radiative-convective interface. Thus z G [0, h] represents the "radiative zone" 



-4- 



while z G [/i, 1] represents the "convection zone". From here on, h = 0.7. Figure [T] illustrates 
the geometry of the Cartesian system. 




z=l 



z=h 



z=0 



y=Jt/2 



Fig. 1. — Cartesian model geometry and intended correspondence with the spherical case. 
The shaded area marks the convective region, where forcing is applied. The y = and 
y = 7r/2 lines mark the "poles" and the "equator". The system is assumed to be periodic 
with period 7r in the y— direction. 



2.2. Model equations 

For this simple Cartesian approach, we work with the Boussinesq approximation (this 
assumption is dropped in the spherical model but doesn't affect the predicted scalings). 
The background state is assumed to be stratified, steady, and in hydrostatic equilibrium. 
The background density and temperature profiles are denoted by T(z) and p(z) respectively. 
Density and temperature perturbations to this background state (p and T) are then assumed 
to be linearly related: p = —aT{z)T, where olt{z) is the coefficient of thermal expansion. 
For simplicity, we will assume that the background temperature gradient T z is constant 
throughout the domain z G [0, 1], and treat the convection zone as a region where ax — > 
while the radiative zone has ^ (see below). The alternative option of using a constant a-p 
and a varying T z yields qualitatively equivalent scalings for the meridional flows (a statement 
which is verified in Section [3]) although the algebra is trickier. 

The set of equations governing the system in this approximation are the momentum, 
mass and thermal energy conservation equations respectively: 

— + u • Vu + 2e 2 x u = -Vp + Ro 2 (z)Te z + E^V 2 u , 
at 

V-u = , 



where u = (u,v,z) is the velocity field. In these equations, temperature perturbations 
have been normalized to the background temperature difference across the box R Q T Z . The 
governing non-dimensional parameters are: 



Ro(z) = N(z)/fl Q the Rossby number , 



Ej, = — 2 the Ekman number , 



Pr = — the Prandtl number , 



(2) 



where the dimensional quantity N 2 (z) = ax{z)T z g is the square of the Brunt- Vaisala fre- 
quency (g is the magnitude of gravity and is assumed to be constant). The microscopic 
diffusion coefficients v (the viscosity) and k t (the thermal conductivity) are both assumed 
to be constant. In the Sun, near the radiative-convective interface, E u ~ 2 x 10~ 15 and Pr 
~2x KT 6 (Gough, 2007). 

As mentioned earlier, the transition between the model radiative zone and convection 
zone is measured by the behavior of o:t{z), which goes from for z > h to a finite value for 
z < h. This can be modeled in a non-dimensional way through the Rossby number, which 
we assume is of the form: 



where Ro rz is constant. The lengthscale A may be thought of as the thickness of the "over- 
shoot" region near the base of the convection zone, but in practise is mostly used to ensure 
continuity and smoothness of the background state through the tanh function. 

In what follows, we further restrict our study to an "axially" symmetric (d/dx = 0), 
steady-state (d/dt = 0) problem. Within the radiative zone, the nonlinear terms u • Vu 
and u • VT are assumed to be negligible. Within the convection zone on the other hand, 
anisotropic turbulent stresses are thought to drive the observed differential rotation^. We 
model this effect in the simplest possible way, replacing the divergence of the stresses by a 
linear relaxation towards the observed convection zone profile: 




(3) 



u — u. 



■cz 



(4) 



u ■ Vu 



r 



1 Note that for slowly rotating stars, such as the Sun, the direct generation of meridional flows by 
anisotropic stresses is a much weaker effect in the bulk of the convection zone. We neglect it here. 
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where u cz = u cz (y,z)e x and the function u cz (y,z) models the observed azimuthal velocity 
profile in the solar convection zone. This is analogous to the prescription used by Spiegel 
& Bretherton (1968) in their study of the effect of a convection zone on solar spin-down, 
although in their model the convection zone was not differentially rotating. The dimension- 
less relaxation timescale r can be thought of, for example, as being of the same order of 
magnitude as the convective turnover time divided by the rotation period. It is modeled as: 



-\z) = 



A 



1 + tanh 



h 



A 



(5) 



Note that in the real solar convection zone, r varies by orders of magnitudes between the 
surface (r ~ 10 -3 ) and the bottom of the convection zone (r ~ 1). Here, we assume that A 
is constant for simplicity. 



We adopt the following profile for u cz (y, z): 



u cz (y,z) 



Un(z) 



1 + tanh 



h 



iky 



u cz (z)e 



iky 



U (z) = U (h) + S(z-h) , 



(6) 



where k = 2 to match the equatorial symmetry of the observed solar rotation profile. The 
tanh function once again is merely added to guarantee continuity of the forcing across the 
overshoot layer. The function Uq(z) describes the imposed "vertical shear", and is for sim- 
plicity taken to be a linear function of z. If Uo{h) = 0, the forcing effectively vanishes at 
the base of the convection zone. If Uq{K) ^ on the other hand, a strong azimuthal shear 
is forced at the interface. The observed solar rotation profile appears to be consistent with 
U Q {h) and S both being non-zero (and of the order of 0.1, although since we are studying a 
linear problem, the amplitude of the forcing is somewhat irrelevant). Note that if S = the 
forcing velocity u cz has zero vorticity. 

Finally, the observed asphericity in the temperature profile is negligible in the solar 
convection zone; this is attributed to the fact that the turbulent convection very efficiently 
mixes heat both vertically and horizontally. We model this effect as: 



u-VT^ -D(z)V 2 T , 

where the turbulent heat diffusion coefficient is modeled as 

' z — h 



v ; 2 



1 + tanh 



(7) 



(8) 



and thus vanishes beneath the overshoot layer. We will assume that the diffusion timescale 
1/Dq (in non-dimensional units) is much smaller than any other typical timescale in the 
system (D 3> 1). 
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Projecting the remaining equations into the Cartesian coordinate system, and seeking 
solutions in the form q(y, z) = q(z)e tky for each of the unknown quantities yields 



Note that as required, the imposed forcing term drags the fluid in the azimuthal direction: 
for t — > 0, u — > u cz in the convection zone. The meridional flows v and w on the other 
hand are generated by the y— component of the Coriolis force and by mass conservation 
respectively (the essence of gyroscopic pumping, see Mclntyre 2007). 

We now proceed to solve these equations to gain a better understanding of the meridional 
flows and their degree of penetration into the radiative zone below. We use a dual approach, 
solving these equations first analytically under various limits, and then exactly using a simple 
Newton-Raphson-Kantorovich (NRK) two-point boundary value algorithm. The analytical 
approximations yield predictions for the relevant scalings of the solutions in terms of the 
governing parameters (an in particular, the Ekman number and the Rossby number) which 
are then confirmed by the exact numerical solutions. 



Although this limit is not a priori relevant to the physics of the solar interior, we begin 
by studying the case of an unstratified region, setting Ro rz = (in this case, the thermal 
energy equation can be discarded). This simpler problem, as we shall demonstrate, contains 
the essence of the problem. 

In order to find analytical approximations to the solutions, we solve the governing 
equations separately in the convective zone and in the radiative zone. At this point, it 
may be worth pointing out that in the unstratified case, the nomenclatures "convective" 
and "radiative" merely refer to regions which respectively are and are not subject to the 
additional forcing. 




(9) 



2.3. 



The unstratified case 
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We assume that the transition region is very thhfl. In this case, r 1 = A for z > h 
while t -1 = for z < h. Similarly, u cz (z) = Uo(h) + S(z — h) in the convection zone 
while u cz {z) = in the radiative zone. Once obtained, the solutions are patched at the 
radiative-convective interface. 



2.3.1. Solution in the convection zone 

In the convection zone, the equations reduce to 

-2v = -A(u - u cz ) , 
2u = —ikp — Av , 

= -f-Aw, 
dz 

ikv + ^ = 0, (10) 
dz 

where we have neglected the viscous dissipation terms in favor of the forcing terms since 
E„ C A for all reasonable solar parameters. Combining them yields 

dz* ~ 4 + A* W + 2t 4 + A2 dz ■ {U} 

This second-order ordinary differential equation^! for w(z) suggests the introduction of a new 
lengthscale 

6 = ^A—> (12) 

so that the general solution to tTTUj) is 

2iS 



w{z\ 



Ae z/5 + Be~ z/S 



kA ' 

u(z) = u cz (z) - ^ [Ae^ s - Be^ s ] , 
p(z) = -^-u cz {z) - 5A [Ae z ' s - Be~ z ' 5 ] . (13) 

IK 



2 More precisely, A <C Ej/ 2 , see Section [2.3.51 

3 The original order of the system is much reduced in the convection zone since we ignored the effect of 
viscous terms there. 
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The constants A and B are integration constants which must be determined by apply- 
ing boundary conditions (at z — 1) and matching conditions (at z = h). Note from the 
■u-equation that the actual rotation profile approaches the imposed (observed) profile u cz 
provided A and B tend to 0, or when A ^> 2 (in which case 5 — > 1/k). 



2.3.2. Solution in the radiative zone with stress-free lower boundary, and matching 
In the radiative region, the equations reduce to 

2u = —ikp , 

ikv + = , (14) 
dz 

if we neglect viscous stresses in both y and z components of the momentum equation. Note 
that viscous stresses in the x— equation cannot be dropped since they are the only force 
balancing the Coriolis force. These equations are easily solved: 

P(z) = Prz , 

„ / \ ik 
u{z) = -yPrz , 

iE v k z 



v(z) = —p 



W(Z) = W rz —PrzZ , (15) 

where p rz and w TZ are two additional integration constants. Here, we recover the standard 
Taylor- Proudman constraint where in the absence diffusion or any other stresses, the velocity 
must be constant along the rotation axis (here, e z ); in the limit E u — > 0, u(z) and w(z) 
become independent of z, while v{z) — > 0. 

We are now able to match the solution in the radiative zone to that of the convection 
zone. The two constants p TZ and w rz form, together with A and B, a set of 4 unknown 
constants which are determined by application of boundary and matching conditions. Since 
we have neglected viscous effects in the convection zone, we cannot require any boundary or 
matching condition on the horizontal fluid motions. On the other hand, we are allowed to 
impose impermeability w = at the surface (z = 1) and at the bottom (z — 0). Moreover, 
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we request the continuity of the radial (vertical) velocity and of the pressure at the interface 
[z — h). Applying these conditions yields the set of equations 



Ae h/S + Be~ h ' 5 - 



2iS 
~kA 



W r 



W T7 

E u t 
4 







p Iz h 



- U (h) - 5A [Ae h ' s - Be- h ' 5 ] = p rz , 
Ae i/s + Be -i/s _ = Q 



kA 



(16) 



which have the following solution for A and B: 



A 



¥, u hk s j T , 2iS 
2i U °W + kA 


1 


_ e (l-h)/8 


1 + 


E fh5A~\ 




e h/s [! Ef h 5A} 




[1 + 


E fhSA 





B 



2iS 

~kX ( 



l/S _ Ae 2/S _ 



(17) 



These can be substituted back into ( ]T3l) to obtain the meridional flow velocities in the 
convection zone. While the exact form of A and B are not particularly informative, we 
note that in the limit S = (i.e. the forcing velocity has no azimuthal vorticity), both A 
and B scale as E u . This implies that the amplitude of meridional flows everywhere in the 
solar interior scales like E u (even in the convection zone). The physical interpretation of 
this somewhat surprising limit is discussed in Section I2.3.4[ but turns out to be of academic 
interest only (Section 12.3.51) . 

When S ^ then A and B are of order S/ kA in the convection zone regardless of the 
Ekman number, and respectively tend to 



.4 



B 



2iS l-e^-W* 
~ ~kA e h / s - t^- h )l & 
2iS 1 - e^-V/ 6 

~k~A e -h/& _ e (h-2)/S 



+ 0{E V ) , 
+ 0(E„) , 



18) 



as E u — > 0. This implies that w is of order S/kA in the convection zone. Since significant 
flows are locally generated, one may reasonably expect a fraction of the forced mass flux to 
penetrate into the lower region, especially in this unstratified case. 

Using f TlTl) in ( fl6|) . solving for p rz , then plugging p rz into ( |T5l) . we find that the general 
expression for w(z) in the radiative zone z G [0, h] is 



iE„k 3 



w(z) 



U {h) + 



cosh((l - h)/S) - 1 
sinh((l - h)/5) 



SS)z . 



(19) 
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This implies that only a tiny fraction of the large mass flux circulating in the convection 
zone actually enters the radiative zone. Instead, the system adjusts itself in such a way as 
to ensure that most of the meridional flows return above the base of the convection zone. 

We now compare this a priori counter-intuitiv^l analytical result with exact numerical 
solutions of the governing equations. The numerical solutions were obtained by solving 
(Q for Ro rz = (unstratified case), and are uniformly calculated in the whole domain (i.e. 
there is nothing special about the interfacial point z = h). The boundary conditions used are 
impermeable boundary conditions at the top and bottom of the domain for w, and stress-free 
boundary conditions for u and v. 

In Figure [2] we compare numerical and analytical solutions for w(z), in a case where the 
forcing function parameters are A = 10 -4 , A = 10, Uq{K) = and S = 1, for four values 
of the Ekman number. The analytical solution is described by equations (IT51) and (|17|) (for 
the convection zone) and (fT5l) (for the radiative zone). As E v — > the numerical solution 
approaches the analytically derived one, confirming in particular that w(z) oc zE v in the 
radiative zone. The convection zone solution is also well-approximated in this case by the 
analytical formula. 

A full 2D visualization of the flow for E v = 10 -4 but otherwise the same governing 
parameters is shown in Figure [3j This figure illustrates more clearly the fact that the 
meridional flows are negligible below the interface, and mostly return within the convection 
zone. Note that given our choice of the forcing function u cz (y,z) oc cos(2y), the induced 
Coriolis force does not vanish at y = or y = n (the "poles"). This explains why the 
meridional flows apparently cross the polar axis in this simple model. This is merely a 
geometric effect: in a true spherical geometry the forcing azimuthal velocity u cz (r, 9) would be 
null at the poles, and meridional flows cannot cross the polar axis. More realistic calculations 
in spherical geometry are discussed in Section [3J 

Finally, it is interesting to note that the analytical solution for the azimuthal velocity 
u(z) exhibits a "discontinuity" across the base of the convection zone, which tends to 

- if/.") = (jjjf + ^) (Ae»<* - Be-"") (20) 

as E u — > 0. The numerical solutions of course are continuous, but the continuity is only 
assured by the viscosity in the system (in the y— direction) and the fact that the overshoot 
layer depth is finite. This is shown in Figure HI together with a comparison of the numerical 
solutions with the analytical solution, again confirming the analytical approximation derived. 



4 but a posteriori obvious, see Section f2. 3. 41 



0.0 0.2 0.4 0.6 0.8 1.0 

z 



Fig. 2. — Numerical (solid) and analytical (dashed) solutions for |it)(z)|, in the case of a 
stress-free bottom boundary. From the uppermost to lowermost curves, E v = 10 -3 , 10 -4 , 
10~ 5 and 10~ 6 respectively, confirming the analytical scaling that w(z) oc E u z in the radiative 
zone while w(z) becomes independent of E u in the convection zone. These solutions were 
obtained with forcing defined by the parameters A = 10~ 4 , A = 10, Uq(z) = S(z — h) and 
S = 1. 

This highlights another and equally a priori counter-intuitive property of the system: 
the value of w rz in the radiative zone is markedly different from the imposed u C z{h) = Uo{h) 
at the interface: 

u rz = U (h) + ^ (Ae h ' 5 - Be- h / s ) . (21) 

Hence, even if the imposed differential rotation is exactly at the radiative-convective 
interface (as it is the case in the simulation presented in Figure H] since U (z) = S(z — h)), 
a large-scale latitudinal shear measured by u rz may be present in the radiative zone, as 
illustrated in Figure [3j This shows that the propagation of the azimuthal shear into the 
radiative zone is non-local (i.e. does not rely on the presence of shear at the interface), and 
is instead communicated by the long-range pressure gradient. 
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Fig. 3. — 2D visualization of the flow for E v = 10~ 4 , in the case of a stress-free bottom 
boundary. Shown as solid and dotted line respectively are linearly spaced streamlines of 
counter-clockwise and clockwise meridional flows. As predicted, the flows appear to return 
entirely within the convection zone and carry a negligible mass flux into the radiative zone. 
Meanwhile the azimuthal velocity (it) as displayed in the filled contours is constant along 
the rotation axis (z— axis) below the interface (z — h — 0.7), but is strongly sheared at the 
interface. This solution was obtained with forcing defined by the parameters A = 10~ 4 , 
A = 10, U (z) = S(z - h) and S = 1, as in Figured 
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Fig. 4. — Numerical (solid) and analytical (dashed) solutions for u(z), in the case of a 
stress-free bottom boundary. From the lowermost to uppermost curves, E v = 1CT 3 , 1CT 4 , 
10 -5 and 10~ 6 respectively, confirming that u(z) tends to a constant in the radiative zone, 
while sustaining a finite discontinuity at the radiative-convective interface (z = h = 0.7). 
These solutions were obtained with forcing defined by the parameters A = 10~ 4 , A = 10, 
Uq(z) = S(z — h) and S — 1, as in Figure 



2.3.3. Solution in the case of no- slip bottom boundary 

The stress-free bottom boundary conditions studied in the previous Section are at first 
glance the closest to what one may expect in the real Sun, where the "bottom" boundary 
merely represents the origin of the spherical coordinate system. However, let us now explore 
for completeness (and for further reasons that will be clarified in the next Section) the case 
of no-slip bottom boundary conditions. 

When the lower boundary is a no-slip boundary, the nature of the solution in the whole 
domain changes. This change is induced by the presence of an Ekman boundary layer, 
which forms near z — 0. Just above the boundary layer, in the bulk of the radiative zone, 
the solution described in 12.3.21 remains valid. However, matching the bulk solution with the 
boundary conditions can no longer be done directly; one must first solve for the boundary 
layer dynamics to match the bulk solution with the boundary conditions across the boundary 
layer. This is a standard procedure (summarized in Appendix A for completeness), and leads 
to the well-known "Ekman jump" relationship between the jump in w(z) and the jump in 
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u(z) across the boundary layer: 

Mbuik - u(0) = jE^ 2 (w hnlk - w(0)) . (22) 

By impermeability, w(0) = 0. Moreover, by assuming that the total angular momentum 
of the lower boundary is the same as that of the convection zone, we require that u(0) = 0. 
Meanwhile, wtmik — u rz and uvik = w TZ in the notation of equation (fT5|) . So finally, for 
no-slip boundary conditions, we simply replace the impermeability condition (w rz = 0) in 
(□ED by 

u xz = jE;^ 2 w rz , (23) 

and solve for the unknown constants A, B, w TZ and p rz as before. 

The exact expressions for the resulting integration constants A and B are now slightly 
different from those given in (iT7|) . but are without particular interest. However, it can be 
shown that they have the same limit as in the stress-free case as E u — > (with 5 ^ 0). This 
implies that the meridional flows driven within the convection zone, in the limit E u — > 0, 
and with 5* 7^ 0, are independent of the boundary condition selected at the bottom of the 
radiative zone. However, we now have the following expression for w(z) in the bulk of the 
radiative zone: 

which has one fundamental consequence: the amplitude of the flows allowed to penetrate into 

1 /2 

the radiative zone is now of order Ej instead of being 0(E U ). This particular statement is 
actually true even if S = 0, although in that case both convection zone and radiative zone 

1 /2 

flows scale with EJ . 

Figure [5] shows a comparison between the approximate analytical formula and the nu- 
merical solution for the same simulations as in Figure [2, but now using no-slip bottom 
boundary conditions. For ease of comparison, the results from the stress-free numerical 
simulations (for exactly the same parameters) have also been drawn, highlighting the much 
larger amplitude of the meridional flows down-welling into the radiative zone in the no-slip 
case, and their scaling with Ej . Figure [6] shows an equivalent 2D rendition of the solution, 
and illustrates the presence of large-scale mixing in the bulk of the radiative zone when the 
bottom-boundary is no-slip. 
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Fig. 5. — Numerical (solid) and analytical (dashed) solutions for in the case of a 

no-slip bottom boundary. From the uppermost to lowermost curves (as seen in the radiative 

zone), E v = 10~ 3 , 10 -4 , 10~ 5 and 10~ 6 respectively, confirming the analytical scaling that 

1/2 

w(z) oc Ej in the radiative zone while w(z) becomes independent of E v in the convection 
zone. These solutions were obtained with forcing defined by the parameters A = 10~ 4 , 
A = 10, Uo(z) = S(z — h) and S — 1, as in Figure El For comparison, the previous 
simulations with stress-free bottom boundary, for the same parameters, are shown as dotted 
lines. 

2.3.4- Physical interpretation 

The various sets of solutions derived above can be physically understood in the following 
way. Let us first discuss the solution in the convection zone. In the limit where u cz (y,z) is 
independent of z (equivalently, S = 0), the azimuthal (x—) component of the vorticity of 
the forcing is zero. In that case there is no injection of x— vorticity into the system aside 
from that induced in the viscous boundary layers, and the amplitude of the meridional flows 
generated in the convection zone scales with E u . This limit is somewhat academic in the 
case of the Sun, however given the observed rotation profile (see also Section H2.3. 51) . When 
S 7^ 0, the amplitude of the induced meridional flows in the convection zone scales linearly 
with S and is independent of viscosity. 

In the radiative zone, the Taylor- Proudman constraint enforces invariance of the flow 
velocities along the rotation axis, except in regions where other forces balance the Coriolis 
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Fig. 6. — The same as for Figure [3] but for no-slip boundary conditions. The Ekman layer 
near the lower boundary is clearly visible. For ease of comparison, the same streamlines are 
shown in the two plots. The two figures illustrates how the nature of the lower boundary 
condition influences the mass flux through the radiative zone. 
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force. In the non-magnetic, unstratified situation discussed in the two previous sections, 
the only agent capable of breaking the Taylor-Proudman constraint are viscous stresses, 
which are only significant in two thin boundary layers: one right below the convection zone 
and the other one near the bottom boundary. These two layers are the only regions where 
flows down-welling into the radiative zone are allowed to return to the convection zone. The 
question then remains of what fraction of the mass flux entering the radiative zone returns 
within the upper Ekman layer, and what fraction returns within the lower Ekman layer. The 
latter, of course, permits large-scale mixing within the radiative zone. 

In the first case studied, the bottom boundary was chosen to be stress-free. This nat- 
urally suppresses the lower viscous boundary layer so that the only place where flows are 
allowed to return is at the radiative-convective interface. As a result, only a tiny fraction of 
the mass flux penetrates below z = h, and the turnover time of the remaining flows within 
the radiative zone is limited to a viscous timescale of the order of 1/E U Q & . 

Following this reasoning, we expect and indeed find quite a different behavior when the 
bottom boundary is no-slip. In that case viscous stresses within the lower boundary layer 

1 /2 

break the Taylor-Proudman constraint and allow a non-zero mass flux (of order ) to 
return near z — 0. This flow then mixes the entire radiative zone as well, with an overall 

1 /2 

turnover time of order of 1/~EJ f2 Q (in dimensional units). 

To summarize, in this unstratified steady-state situation, the amount of mixing induced 
within the radiative zone by convective zone flows depends (of course) on the amplitude of 
the convection zone forcing, but also on the existence of a mechanism to break the Taylor- 
Proudman constraint somewhere within the radiative zone. That mechanism is needed in 
order to allow down-welling flows to return to the convection zone. But more crucially, this 
phenomenon implies that the dynamics of the lower boundary layer entirely control the mass 
flux through the system. 

Here, we studied the case of viscous stresses only. One can rightfully argue that there 
are no expected "solid" boundaries in a stellar interior and that the overall behavior of the 
system should be closer to the one discussed in the stress-free case than the no-slip case. 
However, we chose here to study viscous stresses simply because they are the easiest available 
example. In real stars viscous stresses are likely to be negligible compared with a variety 
of other possible stresses: turbulent stresses at the interface with another convection zone, 
magnetic forces, etc. Nevertheless, these stresses will play a similar role in allowing flows to 
mix the radiative zone if they become comparable in amplitude with the Coriolis force, and 
help break the Taylor-Proudman constraint. This issue is discussed in more detail in Section 
H 
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2.3.5. The thickness of the overshoot layer 

Before moving on to the more realistic stratified case, note that this unstratified system 
holds one final subtlety. In all simulations presented earlier, the overshoot layer depth was 
selected to be very small - and in particular, smaller than the Ekman layer thickness. In 
that case, the transition in the forcing at the base of the convection zone is indeed close to 
being a discontinuity, and the analytical solutions presented in Sections 12.3.11 and 12.3.21 are 
a good fit to the true numerical solution. 

In the Sun, the overshoot layer depth A is arguably always thicker than an Ekman 

1 12 

lengthscale ~E U . When this happens, the solution "knows" about the exact shape of the 
forcing function within the transition region, and therefore depends on it. This limit turns 
out to be rather difficult to study analytically, and since in the case of the Sun we do not 
know the actual profile of t~ 1 (z), there is little point in exercise anyway. 

We can explore the behavior of the system numerically, however, for the profile t~ 1 (z) 
discussed in equation (jSj), in the limit where A > . The example for which this effect 
matters the most is the somewhat academic limit where S = in the convection zone, but 
Uo(h) 0. In this case, the asymptotic analysis predicts that the meridional flow amplitudes 

1 /2 

are 0(E U ) in both the convection zone and in the radiative zone for the no-slip case. We 

1 /2 

see in Figure [7] that this is indeed the case in simulations where A <C Fj u ■ However, when 
the overshoot thickness is progressively increased and becomes larger than the Ekman layer 
thickness, the amplitude of the meridional circulation in the convection zone is no longer 

1 /2 

0(E„ ) but much larger. Meanwhile, the scaling of the radiative zone solutions with Fj u 
remain qualitatively correct. The difference with the analytical solution in the convection 
zone can simply be attributed to the fact that when the system knows about the shear within 
the overshoot layer the limit S = is no longer relevant. 



2.4. The stratified case 

While the previous section provides interesting insight into the problem, notably on the 
role of the Taylor-Proudman constraint, we now move to the more realistic situation where 
stratification plays a role in the flow dynamics. In this section, we generalize our Cartesian 
study to take into account the stratification of the lower region (Ro rz ^ 0). For this purpose, 
we go back to studying the full system of equations (jHJ). As before, we first find approximate 
analytical solutions to derive the overall scaling of the solutions with governing parameters, 
and then compare them to the full numerical solutions of (19"1). The analytical solutions are 
obtained by solving the system in the convective region and radiative region separately, and 
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Fig. 7. — Comparison of numerical simulations (solid lines) with analytical prediction 
(dashed line) for a no-slip bottom boundary, for E u = 10~ 6 , and forcing functions defined 
with A = 10, S — and U (h) = 1. The three numerical solutions are obtained for various 
values of the overshoot layer depth: from lowermost to uppermost curves (as seen in the con- 
vective zone), A = lOE^ 2 , E l J 2 , and 0.1E, 1 / 2 . The analytical solution assumes an infinitely 

thin overshoot layer and is therefore independent of A. Note that the analytical solution 

1 li 

in the convection zone is only a good approximation to the true solution if A EJ . The 
overall scalings in the radiative zone, however, are preserved. 

matching them at z — h. 



2.4-1- Convection zone solution 

The equations in the convection zone are now given by 

-2v = -A(u - u cz ) , 
2u = —ikp — Av , 

dz 
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ikv + ^ = , (25) 
dz 

where we have assumed that D ^> E^/Pr. Eliminating variables one by one yields the same 
equation for w(z) in the convection zone as before fTTTT) . as well as a Poisson equation for T 
once w is known. The solutions are then as ({TBI) , together with 

T i ,)=T«P + T 1 e-+ + L__ J + _, (26) 

where the integration constants T and remain to be determined. For the sake of analytical 
simplicity, we will assume that D Q ^> 1 in all that follows (i.e. very large thermal diffusivity 
in the convection zone), and thus neglect the third and fourth terms in ( 1261) . This limit is 
relevant for the Sun. 



2.4-2. Radiative zone solution 

The radiative zone equations are now 

/d 2 u 

- 2v = E v ( — - k z u 



dz 2 
2u = —ikp 

= -g+Ro£,(z)T 

. E„ ( d 2 f 



ikv + ^ = , (27) 
dz 



and can be combined to yield 



d A u ,,/ PrRo r 2 \ d 2 u , , PrRo 



1 + ^ + fc*l^p-fl = 0, (28) 



dz 4 V 4 / d -2 2 

and similarly for T. The characteristic polynomial is 



with solutions 



(\ 2 -k 2 )(\ 2 -^^k 2 )=0, (29) 



± Ai = ±k , 

±A 2 = ±VP^^k . (30) 
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These solutions are the same as those presented in Paper I, and will be referred to as the 
"global-scale" mode and the thermo-viscous mode respectively. Note that here, A 2 corre- 
sponds to k 2 in Paper I. 

In this steady-state study, the quantity A2 summarizes the effect of stratification. It is 
important to note that it contains information about the rotation rate of the star as well as 
the Prandtl number, in addition to the buoyancy frequency. If A2 <C 1, then the thermo- 
viscous mode essentially spans the whole domain: the system appears to be "unstratified" , 
and is again dominated by the Taylor- Proudman constraint. On the other hand, if A2 ^> 1 
then the flows only penetrate into the radiative zone within a small thermo-viscous boundary 
layer of thickness I/A2 as a result of the strong stratification of the system. The Taylor- 
Proudman constraint is irrelevant in this limit, since the magnitude of the buoyancy force is 
much larger than that of the Coriolis force. 

The calculation above was made in the limit where the viscous terms in the latitudi- 
nal and radial components of the momentum equation are discarded. Paper I shows that 
two additional Ekman modes are also present if they are instead kept. By analogy with 
the unstratified case, we expect that these Ekman modes do not influence the solution for 
stress-free boundary conditions, but that additional care must be taken for no-slip boundary 
conditions. 

Note that the equation for w instead simplifies to 

d 2 w 2 PrRo rz 2 A 

d^ = k —^ W > (31) 

and similarly for v (i.e. both equations are only second order in z, and only contain the 
thermo-viscous mode). The radiative zone (z G [0, h]) solutions are now 

u(z) = u x e kz + u 2 e~ kz + u 3 e X2Z + u A e~ Mz , 
v(z) = ^-{k 2 -\l)[u^ z + u,e-^ z ] , 

(2) = —ik^- ^ 2 ~ — [u 3 e X2Z - u 4 e- X2Z ] , 
2 A 2 

2 

p(z) = — - \uie kz + u 2 e~ kz + u 3 e X2Z + u 4 e~ X2Z ] , 
ik L J 

f{z) = — =- \k Ul e kz - ku 2 e~ kz + \ 2 u 3 e x * z - X 2 u 4 e~ X2Z ] . (32) 

ikRo rz 

where the 4 constants {ui} i= i^ are integration constants, to be determined. 
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2.4-3. The stratified stress-free case 

We now proceed to match the solutions in the two regions, assuming stress-free boundary 
conditions near the lower boundary. Since there are in total 8 unknown constants (including 
A, B, T and T\ from the convection zone solution and {ui\i = \^ from the radiative zone 
solution), we need a total of 8 matching and boundary conditions. 

At the lower boundary (z = 0) we take w = du/dz = 0; this condition in turn implies 
that T = 0. At the surface [z = 1), we take as before w = 0, and select in addition T = 0. 
We then need 4 matching conditions across the interface: these are given by the continuity 
of w, p, T and dT/dz. Note that it is important to resist the temptation of requiring the 
continuity of v, since viscous stresses have been neglected in the analytical treatment of the 
y— component of the momentum equation in both radiative and convective zones. Moreover, 
we know that in the unstratified limit, u actually becomes discontinuous at the interface in 
the limit E„ — > 0. Since we expect the stratified solution to tend to the unstratified one 
uniformly as Ro rz — > 0, we cannot require the continuity of u at the interfac^]. 

The equations and resulting solutions for the integration constants are fairly compli- 
cated. The most important ones are reported in the Appendix B for completeness, and are 
used to justify mathematically the following statements: 

• In the limit of Ro rz — > 0, we find as expected that the solutions uniformly tend to 
the unstratified solution summarized in equations f[T3"j) . (JT5]1 . and (fT7|) . Indeed, in 
that case A2 — > and the thermo-viscous solution spans the whole radiative zone 
(mathematically, it tends to the linear solution found in the unstratified case). 

• In the strongly stratified case (defined as A2 k), as described earlier, w in the radia- 
tive zone decays exponentially with depth on a lengthscale I/A2, with an amplitude 
which scales as E u /\2- The flows are therefore very strongly suppressed, and return to 
the convection zone within a small thermo-viscous layer. Note that E„/A2 = Ra -1 ^ 2 
where Ra is the usually defined Rayleigh number. 

The two limits are illustrated in Figure [BJ which shows the numerical solution to (Q for 
two values of the Rossby number Ro rz , but otherwise identical parameters. In the strongly 
stratified limit (A 2 = 10, using Pr = 0.01 and Ro rz = 10 2 ) we see that the solution decays 
exponentially below the interface, with an amplitude which scales as Fi u /\ 2 as predicted 
analytically. In the weakly stratified case (A2 = 0.1, using Pr = 0.01 and Ro rz = 1) the 
solution tends to the unstratified limit and scales as E„z. 



5 a fact which is again only obvious in hindsight 
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Fig. 8. — Numerical solutions of ([9]) with the following parameters: A = 0.01, A = 10, 
S = 0, U = 1, Pr = 0.01 and D = 10. Stress-free bottom boundary conditions are used. 
The solid lines correspond to the "strongly" stratified case with Ro rz = 10, with E u = 10~ 5 
and 10 -6 for the top and bottom curves respectively. The dashed lines correspond to the 
"weakly" stratified case, with Ro rz = 1, with E u = 10~ 5 and 10~ 6 for the top and bottom 
curves respectively. Note that for k = 2, A2 is simply equal to Pr 1 / 2 Ro rz . For comparison, 
the unstratified case (Ro rz = 0) is shown as dotted lines. At these parameters and with these 
boundary conditions, Ro rz = 1 already belongs to the weakly stratified limit. 

2.4-4- Matching in the no-slip case 

By analogy with the previous section, we expect to recover the unstratified limit when 

1/2 

A2 — > 0, so that w(z) oc FjJ in this no-slip case. In the strongly stratified limit on the 
other hand, the amplitude of the flows decays exponentially with depth below the interface 
as a result of the thermo-viscous mode and is negligible by the time they reach the lower 
boundary. In that case, we do not expect the applied lower boundary conditions to affect the 
solution, so that the scalings found in the strongly stratified limit with stress-free boundary 
conditions should still apply: w(z) oc E u /\ 2 . 

These statements are verified in Figure M There, we show the results of a series of 
numerical experiments for no-slip boundary conditions where we extracted from the simu- 
lations the power a in the expression «; oc E", and plotted it as a function of stratification 
(A 2 ). To do this, we integrated the solutions to equations (jSJ) for the following parameters: 
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A = 0.01, A = 
4 values of E,, 



10, 5 = 1, U (h) = 1, Pr = 0.01 and D = 10 and calculated w{z = 0.5) for 
10~ 6 , 10~ 7 , 10~ 8 and 10~ 9 . We estimated a by calculating the quantity 



a 



log 



W(Z 



0.5, E, 



io- 



10 



w (z 



0.5, E u = 1Q- 7 ) 



(33) 



for the (Ej, 



10 



-6 



E u = 10 7 ) pair (diamond symbols) and similarly for the pairs (E„ 



10~ 7 , E u = 10 -8 ) (triangular symbols) and (E„ = 10~ 8 , E u = 10~ 9 ) (star symbols). In the 
weakly stratified limit (A2 — > 0), we find that a — ^ 1/2 while in the strongly stratified limit 
(A2 S> 1), a — > 1, thus confirming our analysis. The transition between the two regimes 
appears to occur for slightly lower-than expected values of A 2 , namely 0.1 instead of 1. 
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Fig. 9. — This figure shows the power a in the expression w cx E", as a function of A2 (see 
main text for detail). In the weakly stratified limit, a — > 1/2 while in the strongly stratified 
limit a — > 1 as predicted analytically. This calculation was done for no-slip boundary 
conditions, and the following parameters were held constant: A = 0.01, A = 10, S — 1, 
U (h) = 1, Pr = 0.01 and D = 10 

A final summary of our findings for the stratified case together with its implications for 
mixing between the solar convection zone and the radiative interior, is deferred to Section 
HI There, we also discuss the consequences in terms of mixing in other stars. But first, we 
complete the study by releasing some of the simplifying assumptions made, and moving to 
more realistic numerical solutions to confirm our simple Cartesian analysis. 
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3. A "solar" model 



In this section we improve on the Cartesian analysis by moving to a spherical radiative- 
convective model. The calculations are thus performed in an axisymmetric spherical shell, 
with the outer radius r out selected to be near the solar surface, and the inner radius r in 
somewhere within the radiative interior. This enables us to gain a better understanding of 
the effects of the geometry of the system on the spatial structure of the flows generated. In 
addition, we use more realistic input physics in particular in terms of the background stratifi- 
cation, and no longer use the Boussinesq approximation for the equation of state. We expect 
that the overall scalings derived in Section [2] still adequately describe the flow amplitudes in 
this new calculation. However, the use of a more realistic background stratification adds an 
additional complication to the problem: the background temperature/density gradients are 
no longer constant, so that the measure of stratification A2 varies with radius (see Figure fTO} 
for an estimate of A2 in the Sun). This aim of this section is therefore to study the impact 




Fig. 10. — Variation of A2/A; =Pr°' 5 A^/f2 in the Sun, as determined from Model S of 
Christensen-Dalsgaard et al. (1996). The Prandtl number Pr is calculated using Model 
S together with the formulae provided by Gough (2007) for the microscopic values of the 
viscosity v and the thermal conductivity kt (see also Garaud & Garaud, 2008). The inset 
zooms into the region near the base of the convection zone, which is the only region of the 
radiative zone where A2 < 1 (aside from r — > 0). 

of both geometry and non-uniform stratification on the system dynamics. 
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3.1. Description of the model 

The spherical model used is analogous to the radiative-zone-only model presented in 
Paper I and described in detail (including the magnetic case) by Garaud & Garaud (2008). 
The salient points are repeated here for completeness, together with the added modifications 
made to include the "convective" region. 

We consider a spherical coordinate system (r, 8, 0) where the polar axis is aligned with 
rotation axis of the Sun. The background state is assumed to be spherically symmetric and 
in hydrostatic equilibrium. The background thermodynamical quantities such as density, 
pressure, temperature and entropy are denoted with bars (as p(r), p(r), T(r) and s(r) 
respectively), and extracted from the standard solar model of Christensen-Dalsgaard et al. 
(1996). Perturbations to this background induced by the velocity field u = (u r , ue, u$) are 
denoted with tildes. In the frame of reference rotating with angular velocity Qq, in a steady 
state, the linearized perturbation equations become 

V • (pu) = , 

2pQ, Q e z x u = -Vp + pg + /V ■ n - pfl Q - 



r(r) 

oc TN 2 r - ~i 
Ur = V ■ [(fkr + Rln Q D(r))VT 



M + L (34) 
p p T 

where c p is the specific heat at constant pressure, kr(r) = pc p Rt is the thermal conductivity 
in the solar interior, II is the viscous stress tensor (which depends on the background viscosity 
v) and g = —g(r)e r is gravity. Note that this set of equations is given in dimensional form 
here, although the numerical algorithm used further casts them into a non-dimensional 
form. Also note that both diffusion terms (viscous diffusion and heat diffusion) have been 
multiplied by the same factor /. This enables us to vary the effective Ekman number 
Ej,(r) = /E® = fv(r) / 'RqQq while maintaining a solar Prandtl number at every radial 
position. As a result, the quantity A 2 used in the simulations and represented in Figure [TU1 
is the true solar value (except where specifically mentioned). 

As in the Cartesian case, we model the dynamical effect of turbulent convection in the 
convection zone through a relaxation to the observed profile in the momentum equation, and 
a turbulent diffusion in the thermal energy equation. The expressions for the non-dimensional 
quantities r(r) and D{r) are the same as in equations §5§ and ([HD with z replaced by r/R Q , 
and h = 0.713 instead of h = 0.7 (Christensen-Dalsgaard et al. 1996). In what follows, we 
take A to be 0.01 (i.e. the overshoot layer depth is 1% of the solar radius) although the 
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choice of A has little influence on the scalings derived. The rotation profile in the convection 
zone u cz is selected to be 

u cz (r,0) =r sin 9Vt cz (9) , (35) 

where 

^cz(^) = ^cq (l — a-2 cos 2 9 — a 4 cos 4 9) , (36) 

with 

a 2 = 0.17 , a 4 = 0.08 , 

^ = 463 nHz , (37) 

2,7V 

which is a simple approximation to the helioseismically determined profile (Schou et al. 1998; 
Gough, 2007). Here, Q eq is the observed equatorial rotation rate. As in Paper I, we finally 
select fl Q to be 

^ = ^ cq (l-f -f^) , (38) 

to ensure that the system has the same specific angular momentum as that of the imposed 
profile u cz (r,6>). 

The computational domain is a spherical shell with the outer boundary located at 
^out = O.9i? . It is chosen to be well-below the solar surface to avoid complications related 
to the very rapidly changing background in the region r > O.95i? . The position of the lower 
boundary will be varied. 

The upper and lower boundaries are assumed to be impermeable. The upper boundary 
is always stress-free, while the lower boundary is assumed to be either no-slip or stress- free 
depending on the calculation. In the no-slip case, the rotation rate of the excluded core is 
an eigenvalue of the problem, calculated in such a way as to guarantee that the total torque 
applied to the core is zero. Finally, the boundary conditions on temperature are selected in 
such a way as to guarantee that V 2 T = outside of the computational domain, as in Garaud 
& Garaud (2008). We verified that the selection of the temperature boundary conditions 
only has a qualitative influence on the results, and doesn't affect the scalings derived. 

The numerical method of solution is based on the expansion of the governing equations 
onto the spherical coordinate system, followed by their projection onto Chebishev polynomi- 
als T n (cos#), and finally, solution of the resulting ODE system in r using a Newton- Raphson- 
Kantorovich algorithm. The typical solutions shown have 3000 meshpoints and 60-80 Fourier 
modes. For more detail, see Garaud (2001) and Garaud & Garaud (2008). 
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3.2. The weakly stratified case 

We first consider the artificial limit of weak stratification. In the following numerical 
experiment, we use the available solar model background state, but divide the buoyancy 
frequency N by 10 3 everywhere in the computational domain (all other background quantities 
remain unchanged). As a result, the new value of A2 in the domain is artificially reduced 
from the one presented in Figure [TU1 by 10 3 , and is everywhere much smaller than one. The 
position of the lower boundary is arbitrarily chosen to be at r- m = 0.35R Q . 

Two sets of solutions are computed for no-slip lower boundary and for stress-free lower 

boundary. Figure [TT] is equivalent to Figure [5j it displays the radial velocity u r as a function 

of radius near the poles (latitude of 80°) for various values of / - in other words, E v - and 

1/2 

clearly illustrates the scalings of u r oc EJ in the radiative zone for the no-slip case, and 
u r oc E v for the stress-free case. 




1 cr 6 f ] 
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Fig. 11. — Vertical velocity at 80° latitude in units of -R©fi© for an artificially weakly 

stratified simulation (where N was uniformly divided by 10 3 everywhere). The solid lines 

show three simulations for / = 10 11 , / = 10 10 and / = 10 9 (from top to bottom) for the 

no-slip case. These correspond to E v = 2 x 10 -4 , E u = 2 x 10~ 5 and E u = 2 x 10~ 6 at the base 

1 12 

of the convection zone respectively, hence showing how u r oc EJ . The dotted lines show 
simulations with stress-free boundary conditions for the same parameters, showing u r oc E v . 
In this calculation the overshoot depth A was selected to be 0.01i? Q , and A = 10. The value 
of D is irrelevant in this very weakly stratified simulation. 
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Figure [12] illustrates the geometry of the flow in both no-slip and stress-free cases for 
/ = 10 9 (which corresponds to an Ekman number near the radiative-convective interface 
of about 2 x 10~ 6 ). The geometrical pattern of the flows observed within the convection 
zone show a single-cell, with poleward flows near the surface and equatorward flows near 
the bottom of the convection zone. Below the convection zone we note the presence of three 
distinct regions: the polar region, a Stewardson layer region (at the tangent cylinder) and 
an equatorial region. Flows within the equatorial region are weak regardless of the lower 
boundary conditions. In the stress-free case, even in the tangent cylinder the flows tend 
to return mostly within the convection zone. If the lower boundary is no-slip on the other 
hand, flows within the tangent cylinder are stronger, although the effect is not as obvious as 
in the Cartesian case because of the anelastic mass conservation equation used here. 




Fig. 12. — Normalized angular velocity (Q/Q eq ) and streamlines solutions to equations (IMI) 
for an artificially weak stratification (see Figure [TT|) . and for / = 10 9 (corresponding to 
E„ = 2 x 10~ 6 at the base of the convection zone). On the left, we show the solution 
with no-slip lower boundary conditions, and on the right the stress-free solution. Dotted 
lines represent clockwise flows, and solid lines counter-clockwise flows. In this calculation, 
the overshoot layer depth A was selected to be 0.01_R Q , and A = 10. The value of D is 
irrelevant in this very weakly stratified simulation. 
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3.3. The stratified case 

Let us now consider the case of a true solar stratification. Since A2 increases rapidly 
with depth beneath the convection zone (from to about 10 in the case of the Sun), the 
radius at which A2 — 1 (7*1 for short) plays a special role: we expect the dynamics of the 
system to depend on the position of the lower boundary r in compared with 7*1 ~ O.55_R . 

This is indeed observed in the simulations, as shown in Figure [131 If > 7"i, then 
A2 < 1 everywhere in the modeled section of the radiative zone. In this case, the dynamics 
follow the scaling for the unstratified case, and depend on the nature of the lower boundary 

1/2 

[u r oc E u if the lower boundary is stress- free, and u T oc Ej if the lower boundary is no-slip). 
On the other hand, if rj n < ri then the flows are strongly quenched by the stratification 
before they reach the lower boundary. As a result, the radial velocities scale with E^/A2 
regardless of the applied boundary conditions. 

The implications of this final result, namely the importance of the location of the stresses 
involved in breaking the Taylor-Proudman constraint in relation to the radius at which 
A2 — 1, are discussed in Section H~3l 

4. Implications of this work for solar and stellar mixing 

4.1. Context for stellar mixing 

The presence of mixing in stellar radiative zones has long been inferred from remaining 
discrepancies between models-without-mixing and observations (see Pinsonneault 1997 for a 
review). The most commonly used additional mixing source is convective overshoot, whereby 
strong convective plumes travel beyond the radiative-convective interface and cause intense 
but very localised (both in time and space) mixing events (Brummell, Clune & Toomre, 
2002). The typical depth of the layer thus mixed, the "overshoot layer", is assumed to be a 
small fraction of a pressure scaleheight in most stellar models. 

A related phenomenon is wave-induced mixing (Schatzman, 1996). While most of the 
energy of the impact of a convective plume hitting the stably stratified fluid below is con- 
verted into local buoyancy mixing, a fraction goes into the excitation of a spectrum of gravity 
waves, which may then propagate much further into the radiative interior. Where and when 
the waves eventually cause mixing (either through mutual interactions, thermal dissipation 
or by transferring momentum to the large-scale flow) depends on a variety of factors. It has 
recently been argued that the interaction of the gravity waves with the local azimuthal veloc- 
ity field (the differential rotation) would dominate the mixing process (Charbonnel & Talon 
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2005), although this statement is not valid unless the wave-spectrum is near-monochromatic. 
For the typically flatter wave spectra self-consistently generated by convection, wave-induced 
mixing has much more turbulent characteristics (Rogers, MacGregor & Glatzmaier 2008), 
and is again fairly localized below the convection zone. 

Mixing induced by large-scale flows comes in two forms: turbulent mixing resulting 
from instabilities of the large-scale flows, and direct transport by the large-scale flows them- 
selves. The former case is the dominant mechanism in the early stages of stellar evolution 
when the star is undergoing rapid internal angular-momentum "reshuffling" caused by ex- 
ternal angular-momentum extraction (disk-locked and/or jet phase, early magnetic breaking 
phase). In these situations, regions of strong radial angular- velocity shear develop, which 
then become unstable and cause local turbulent mixing of both chemical species and an- 
gular momentum. Studies of these processes were initiated by the work of Endal & Sofia 
(1978). Later, Chaboyer & Zahn (1992) refined the analysis to consider the effect of stratifi- 
cation on the turbulence, and showed how this can differentially affect chemical mixing and 
angular-momentum mixing. Finally, Zahn (1992) proposed the first formalism which com- 
bines mixing by large-scale flows and mixing by (two-dimensional) turbulence. In addition 
to the flows driven by angular-momentum redistribution, he also considered flows driven by 
the local baroclinicity of the rotating star, and showed that their effect can be represented 
as a hyper diffusion term in the angular velocity evolution equation. In a quasi-steady, uni- 
formly rotating limit, the flows described are akin to a local Eddington-Sweet circulation. 
His formalism is used today in stellar evolution models with rotation (Maeder & Meynet, 
2000). 

In all cases described above, mixing is either localised near the base of the convection 
zone (overshoot, gravity-wave mixing), or significant only in very rapidly rotating stars (local 
Eddington-Sweet circulations) or stars which are undergoing major angular-momentum re- 
distribution (during phases of gravitational contraction, spin-down, mass loss, etc.). In this 
paper, we have identified another potential cause of mixing, where the original energy source 
is the differential rotation in the stellar convective region: gyroscopic pumping (induced by 
the Coriolis force associated with the differential rotation, see Mclntyre 2007) drives large- 
scale meridional flows which may - under the right circumstances - penetrate the radiative 
region and cause a global circulation. 

This source of mixing is intrinsically non-local to the radiative zone. The simplest way 
of seeing this is to consider a thought-experiment where the radiative-convective interface is 
impermeable: as shown in Paper I, the amplitude of the meridional flows generated locally 
(i.e. below the interface) is then much smaller than the one calculated here. The origin 
of the flows is also clearly independent of the baroclinicity (since the same phenomenon is 
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observed in the unstratified limit), although the flows themselves can be influenced by the 
stratification. This implies that they are not related to Eddington- Sweet flows. Finally, 
contrary to some of the other mixing sources listed above, the one described here does not 
rely on the system being out-of-equilibrium: it is an inherently quasi-steady phenomenon, 
implying that this process is an ideal candidate for "deep mixing" for stars on the Main 
Sequence. 

The process together with the conditions under which strong mixing might occur, are 
now summarized and discussed. 



4.2. Qualitative summary of our results 

The differential rotation observed in stellar convective envelopes (e.g. Barnes et al. 
2005) is thought to be maintained by anisotropic Reynolds stresses, arising from rotationally 
constrained convective eddies (Kippenhahn, 1963). The details of this particular process 
are beyond the scope of this study, but are the subject of current investigation by others 
(Kitchatinov & Riidiger, 1993, 2005; Rempel, 2005). Instead, we have assumed here the 
simplest possible type of forcing which mimics the effect of convective Reynolds stresses in 
driving the system towards a differentially rotating state. Using this model, we then derive 
the expected mixing caused by large-scale meridional flowg^l 

We first found that large-scale flows are indeed self-consistently driven by gyroscopic 
pumping in the convection zone, as expected (Mclntyre, 2007). The amplitude of these flows 
within the convection zone scales roughly as 

V cz ~ rR^AQ) , (39) 

where (Afl) is the observed equator-to-pole differential rotation, R* is the stellar radius, and 
t is, as discussed in Section [2T2l related to the ratio of the convective turnover time divided 
by the rotation period. Note that for the Sun, with (Afl) ~ O.1O , the typical amplitude 
of the corresponding meridional flows would be of the order of 200 r m/s - which doesn't 
seem too unreasonable given the observations of subsurface flows (Giles et al. 1997) and the 
typical values of r in the solar convection zone (see Section I2~2l) . 

Next, we studied how much mixing these flows might induce in the underlying radiative 
zone. In this quasi-steady formalism, we found that the magnitude of convection-zone- 



6 It is worth noting here that while we expect the details of the flow structure and amplitude to be 
different when a more realistic forcing mechanism is taken into account, the overall scalings derived should 
not be affected. 
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driven flows decays exponentially with depth below the radiative-convective interface on the 
lengthscale l 2 , where l 2 = -R /Pr 1//2 Ro rz , as determined in Section [2.4.21 (see also Gilman & 
Miesch, 2004 and Garaud & Brummell, 2008). This penetration corresponds (in the linear 
regime) to a so-called "thermo-viscous" mode. The limit 1% <C R Q corresponds to a strongly 
stratified limit, where the flow velocities are rapidly quenched beneath the convection zone. 
The limit h ^> Rq corresponds to the weakly stratified case, where the thermo-viscous mode 
spans the whole radiative interior and the stratification has little effect on the flow. It is 
important to note that "weakly stratified" regions in this context can either correspond to 
regions with weak temperature stratification (small N), or in rapid rotation, or with small 
Prandtl number - this distinction will be used later. 

The amplitude of the flows upon entering the radiative zone V TZ , together with l 2 , 
uniquely define the global circulation timescale in the interior (roughly speaking, ^/Kz)- 
In the weakly stratified/rapidly rotating limit, we find that the fraction of the meridional 
mass flux pumped in the convection zone which is allowed to enter the radiative zone is 
strongly constrained by Taylor-Proudman's theorem. This theorem, which holds when the 
pressure gradient^ and the Coriolis force are the two dominant forces and are therefore in 
balance, enforces the invariance of all components of the flow velocities along the rotation 
axis. Hence, flows which enter the radiative zone cannot return to the convection zone unless 
the Taylor-Proudman constraint is broken. However, additional stresses (such as Reynolds 
stresses, viscous stresses, magnetic stresses) are needed to break this constraint. As a result, 
two regimes may exist. If the (weakly stratified/rapidly rotating) radiative zone is in pure 
Taylor-Proudman balance, then the system adjusts itself, by adjusting the pressure field, in 
such a way as to ensure that the convection zone flows remain entirely within the convection 
zone. On the other hand, if there are other sources of stresses somewhere in the radiative 
zone to break the Taylor-Proudman balance, then significant large-scale mixing is possible 
since flows entering the radiative zone are allowed to return to the convection zone. Fur- 
thermore, the resulting meridional mass flux in the radiative zone depends rather sensitively 
on the nature of the mechanism which breaks the Taylor-Proudman constraint (see next 
section) . 

The strongly stratified/slowly rotating limit exhibits a very different behavior. Because 
of the strong buoyancy force, the Taylor-Proudman balance becomes irrelevant, the flows 
are exponentially suppressed, and the induced radiative zone mixing is independent of the 
lower boundary conditions. However, note that since N tends to at a radiative-convective 
interface, there will always be a "weakly stratified" region in the vicinity of any convective 



more precisely, the perturbation to the pressure gradient around hydrostatic equilibrium 
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zone. In that region the dynamics described in the previous paragraph apply. 

4.3. Applications to the Sun and other stars 

In the illustrative model studied here, the only stresses available to break the Taylor- 
Proudman constraint are viscous stresses, which are only significant within the thin Ekman 
layer located near an artificial impermeable inner boundary. We do not advocate that this 
is a particularly relevant mechanism for the Sun! However, it is a useful example of the 
sensitive dependence of the global circulation mass flux on the mechanism responsible for 
breaking the Taylor-Proudman constraint. 

In the limit of weak stratification, we found that if the inner boundary is a stress-free 
boundary then the global turnover time within the radiative zone is the viscous timescale. 
This is because stress-free boundary conditions effectively suppress the Ekman layer. On the 
other hand if the boundary layer is no-slip, then the global mass flux through the radiative 
zone is equal to the mass flux allowed to return through the Ekman layer. In that case, and 
according to well-known Ekman layer dynamics, the overall turnover time within the bulk 
of the domain is the geometric mean of the viscous timescale and the rotation timescale 

1 /2 

(1/Ej/ fi ), which correspond to a few million years only. 

Going beyond simple Ekman dynamics, a much more plausible related scenario for the 
solar interior was studied by Gough & Mclntyre (1998). They considered the same mech- 
anism for the generation of large-scale flows within the convection zone, studied how these 
flows down-well into the radiative zone and interact with an embedded large-scale primordial 
magnetic field. They showed that the field can prevent the flows from penetrating too deeply 
into the radiative zone, while the flows confine the field within the interior. In their model, 
this nonlinear interaction occurs in a thin thermo-magnetic diffusion layer, located somewhat 
below the radiative-convective interface. One can therefore see an emerging analogy with 
the dynamics discussed here: in the Gough & Mclntyre model, the field does act as a some- 
what impermeable barrier, and provides an efficient and elegant mechanism for breaking 
the Taylor-Proudman constraint within the radiative zone. The only significant difference is 
that the artificial Ekman layer is replaced by a more convincing thermo-magnetic diffusion 
layer: the mass flux allowed to down-well into the radiative zone, and mix its upper regions, 
is now controlled by a balance between the Coriolis force and magnetic stresses (instead of 
the viscous stresses). With this new balance, they find that the global turnover time for the 
circulation in the region between the base of the convection zone and the thermo-magnetic 
diffusion layer is of the order of a few tens of millions of years (which is still short compared 
with the nuclear evolution timescale). 
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This mixed region is the solar tachocline. By relating their model with observations, 
Gough & Mclntyre were able to identify the position of the magnetic diffusion layer to be 
just at the base of the observed tachocline (around 0.68 ± O.OLR , see Charbonneau et al. 
1999). This turns out to be close enough to the base of the convection zone for the dynamics 
of the radiative region to be weakly-stratified in the sense used in this paper (see Figure [10] 
and Section 13. 3j) so that the meridional flows are indeed able to penetrate, and do so with 
"significant" amplitude (about 10 -5 cm/s) down to the magnetic diffusion layer. However, 
it is rather interesting to note that the Gough & Mclntyre model could not have worked 
had today's tachocline been observed to be much thicker. It is also interesting to note that 
for younger, more rapidly rotating solar-type stars, a much larger region of the radiative 
zone can be considered "weakly stratified" , possibly leading to much deeper mixed regions if 
these stars also host a large-scale primordial field. The implications of these findings for Li 
burning, together with a few other interesting ideas, will be discussed in a future publication. 
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Equation ffl5|) provides the solution "far" from the lower boundary, in the bulk of the 
fluid. Let us refer to the limit of bulk solutions as z — > as &buik(0 + ) (and similarly for 
the other quantities). We now derive the Ekman solution close to the boundary, for the 
unstratified case. Let's study the problem using the stream-function if> with 
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Appendix A: Ekman jump condition 




(1) 



Moreover, let us assume that within 
equations are then approximated by 



the boundary layer, dip/dz 3> kip. The governing 




E, 



dz* ' 
d 3 ^ 



u 



—ikp + E, 



V 



dz 3 
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dp 
d~z 



which simplify to 



d 5 ^ 



-ikE, 



-4 



d 2 ^ 
d7 ' 



(2) 
(3) 



with solutions 



${z) = + VV 3 + ^2e~ X3Z + ^e Mz + ^e~ X4Z 



u(z 



u 



2 



^ie X3Z - ^ 2 e- X3Z 

<^3 <^3 



+ T-^3e 

A 4 



A42 



A. 



— A42 



(4) 
(5) 



where 

A 3 = (l + i)K 1/2 , A 4 = (1 - ^e; 1 / 2 . 

The growing exponentials are ignored to match the solution far from the boundary layer; 
it then becomes clear that u = Wbuik(0 + ), while —ikip = w)buik(0 + ). Requiring no-slip, 
impermeable conditions at z = implies 



u - 



^0 + ^2 + ^4 = 

A 3 -02 + A 4 ?/> 4 = 

= 



-7^2 - T^4 
A3 A 4 



(6) 



which in turn implies 



/ A3 / 
A 4 



^2 



A, 



A, - A 



2 A 3 + A 4 77-1/2 / 

--00 = 2E V 1 ijj . 



u = — — r-^T - ^ = 'M v ^'"ipo • (7) 

&v ^3^4 

This last equation then uniquely relates the limit of the bulk solution w(0 + ) and w(0 + ) as 

2i 



as 



«bulk(0 + 

yielding the standard Ekman jump condition. 



k E; 1 / 2 w hulk (o' 



(8) 



Appendix B: Stratified stress-free solution 

The boundary conditions discussed in Section 12.3.31 imply the following set of equations. 
At z = 0, w = and u z = (alternatively, T = 0): 

= u 3 - u 4 , 
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= kui — ku 2 + A 2 u 3 — X2U4 . (9) 

At z = 1: w = and f = 0: 

= Ae 1 ' 5 + Be- 1 ' 5 - ^ , (10) 

= T e k + T ie ~ k . (11) 

Finally, matching conditions on w, p, T and dT/dz at z — h: 

_ ikE u ^—^-^-u 3 sinh(A 2 /i) = Ae^ + Be""/* - ^ , 
2ui cosh(£;/i) + 2w 3 cosh(A 2 /i) = C/ (/i) + y 5A [Ae^ 5 - Be"^] , 

T e fc/i + T ie - kh = i^- [jfe Ul sinh(Jfe/i) + A 2 w 3 sinh(A 2 /i)] , 

i/cRo rz 

T e fc/l - T ie - M = [jfeV cosh(fc/i) + A 2 u 3 cosh(A 2 /i)l , (12) 

zArRo^, 

where M4 and ii 2 were already eliminated using equations (llip . We now proceed to eliminate 
A, F, To and Ti, which leaves two equations for u\ and u 3 : 

2G [u\ cosh(A;/i) + w 3 cosh(A 2 /i)] — 5hk 2 E v — -w 3 sinh(A 2 /i) , 

2A9 



5S (e {h - 1)/s (l — G) — i) + GU (h) 



A 2 

(F — 1) smh(kh) + ku\ cosh(kh) = — p-u 3 cosh(A 2 /i) — (F — 1) A 2 m 3 sinh(A 2 /i) (13) 

where the functions F(h, k) and h) are geometric factors defined as 

2 



F(/i, jfe) 



]_ _ e 2fe(/i-l) ' 



e (h-2)/6 _ e -h/S 

These equations form a linear system for and w 3 with 
Mi = -Fu 3 , 

55 ( e (^(l_ G )_l) +Gf/o (/,) 

^3 — 75 — 72 , (J-t)) 

2G [-Hcosh(kh) + cosh(A 2 /i)] - 6Ak 2 E v ^f sinh(A 2 /i) 

and where the function H(h, k, A 2 ) is given as 

„ (h , x , A 2 I cosh(A 2 fr) + (F — 1) sinh(A 2 fr) 
1 ' ' A2j fc cosh(fcfr) + (F - 1) sinh(fcfc) • 1 J 
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These rather opaque solutions can be clarified a little by looking at the various relevant 
limits. For weakly stratified fluids A2 — > 0. Then H(h, k, X 2 ) = 0(A|) — > 0, and so 



u 3 



5S (e^-^il - G) - 1) + GU {h) 
2G - 5kk±E u \ 



+ 0(Xl). 



In the limit E,, 



this then becomes 



U3 = - 
3 2 



U (h) - 5S- 



ui = 0(X 2 2 ) , 
cosh((l-A)/5)' 



sinh((l - 

Folding this back into the original solution in the radiative zone then yields 



w(z) 



-ik 3 E u usz 



ik 3 E v 



U (h) - 5S 



1 - cosh((l - h)/S) 
sinh((l - h)/5) 



which is identical to equation (|T9j) . 

In the opposite, strongly stratified limit, A2 3> k. Then we have instead 



H(h,k,X 2 



Xi 



cosh(A2^) 



k 2 cosh(kh) + (F - 1) sinh(fc/i) ' 



so that this time M3 = 0(A 2 ) — > 0, and in the limit E^ — > 



Ml 



2 cosh(fc/i) 



C/ (/i) - 



1 - cosh((l - h)/S) 
sinh((l - h)/S) 



Folding this back into the equation for w(z) in the radiative zone now yields 



w(z 



ik 3 E u sinh(A2^) 
2A2 cosh(A2/i) 



tanh(/c/i) 
tanh(A;(l - h)) 



U (h) - 5S 



1 - cosh((l - h)/S) 
sinh((l - h)/5) 



therefore justifying the scaling discussed in 12.4.31 
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Fig. 13. — Vertical velocity (in units of R & Q Q ) in nine different simulations, at latitude 80°. 
The background stratification in each case is solar, but the position of the lower boundary 
is moved through the radiative zone from 0.65i? Q to O.35_R and 0.1R & . The solid-line plots 
are for no-slip lower boundary conditions while the dotted lines are for stress-free lower 
boundary conditions. Three simulations are shown in each case: (from lowest to highest 
curve) for / = 10 8 , / = 10 9 and / = 10 10 corresponding to E„ = 2 x 1(T 7 to E, = 2 x 1(T 5 . 
The logarithmic scale clearly shows that u r scales with E„ in the radiative zone in the stress- 
free cases for all values of r in while in the no-slip case, r in scales with E u if r in < 0.6, as 
expected from Figure [101 



